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D INTRODUCED 


There are many theories to aid the airfoil design 
process by computational methods, because of the desire to 
reduce the number and cost of wind tunnel tests. A still 
largely unresolved question is the problem of flow sepa- 
ration. Because the classical boundary layer approximation 
cannot be applied to separated flow calculations, engineers 
have tried to overcome this limitation by developing 
viscous/inviscid interaction approaches or to develop direct 
solutions of the Navier-Stokes equation. 

The purpose of this thesis is to demonstrate the 
capability of the viscous/inviscid interaction method by 
applying Cebeci's interactive computer program to a single 
airfoil (FX 63-137) at three low Reynolds numbers and by 
comparing the results with experimental data. 

Chapter II explains the boundary layer theory. The 
boundary layer equations are derived and the turbulence 
models are introduced. Also, this chapter includes the 
prediction of transition boundary layer calculations and 
flow separation. 

Chapter III introduces the interaction methods. 


Three weak interaction methods are explained briefly and the 
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simultaneous method is presented as a strong interaction 


method. 
Chapter IV describes Cebeci's interactive computer 
program. Input/Output data description and JCL files are 


included. Also, the results of the application of Cebeci's 


program are discussed. 


m 


ІІ. BOUNDARY НК ОНО 


А. DERIVATION OF BOUNDARY LAYER EQUATIONS 

Generally, the thickness of the boundary layer increases 
with viscosity, or it is possible to state that it decreases 
with viscosity, or it is possible to state that it qecre вв 
as the Reynolds number increases. From exact solutions of 
the Navier-Stokes equations, it was seen that the 
boundary-layer thickness (d) is proportional to the square 
root of kinematic viscosity (07: 

d(x) о | vx 
U 

where x is the distance from the leading edge of a flat 


plate. Using the local Reynolds number 


Re - U x/y, d(x) w 1 
X Rex 


For simplicity, assume a two-dimensional, steady 





constant - property flow without body forces and leave the 
stresses unspecified so that the results apply to laminar or 
turbulent flow. If simplifications are са. 5е. писце = 
into Navier-Stokes equations, it must be assumed that the 
boundary layer thickness is very small compared with a 
representative linear dimension (L) of the body, ie. d «cm 


In this way, the solutions obtained from the boundary layer 


Ще 


equations are asymptotic and apply to very large Reynolds 
numbers. 


From the Navier-Stokes equations, 








ar ection: a 1 ФР + 1 дбхх + тог MEUM 
OCI, Ox йг! дх ду 
мрест топ: - 4" ор + ] дбуу + 1 дбху зи ду + у ду 
ду (20700 EON 2 x дү 
сі the continuity equation 15: 
ош ко + 0 (233 
2 X ду 
Inserting the "typical" (order-of-magnitude) values, replace 
the dependent variabies as follows; 
Шесто us , Qu ue | ди че 
2y d 2 x уа 1 
Then Eq. (2.1) can be expressed as 
– 1 др + 1 əzx + 1 xy =u 2u*v ðu 
ою (OX P dY 2 x ay 
2 = 
Це бухло Ó xy/^ ue ue (2.4) 
1 1 I 1 1 


where the typical values were written below the terms to 
which they correspond. 

Considering turbulent flow, all stresses are of the 
same order ie. a general stress must be of order 


2. 
P Ue d/l, then Eq. (2.1) becomes 


- 1 др + 1 дбху J ат (4) - u 2210 (72885) 
эш ох г ду 1 дх ду 


2 


where o(d /1) indicates a quantity of the order Gf ааш = 
of 471. 
Considering a laminar flow of Newtonian viscous 


fluid; 





Ó xx = 2U ðu xy =u 2u +2) Ó yy - 24 2v (276! 
OX, ay ðX j ду 


In the Ó xy term, ду/дх is smaller than Qu/oy. 











Therefore, Eq. (2.5) becomes 
2 2 
Нор + )) ди Е + o£ J |- u Qu + v ди (228 
ОУ уу ду“ 1 дх ду 
oimilarly, Eq. (2.2) can be written as: 
- 1 др + 1 30yy «1 2Óxy -u ду + у ду 
сәү че уу P Әх дх ду 
INEO У це р Ued ) ue d ue d. (2.8) 
е фу 197 “үлэ MES p 


If we write all the viscous terms together, 


Әр ud V. ММА 44-4 ue d. (2.9) 


E y 1" uel i 
It is known that (d71Y c^ V/ue1 is laminar flow, so Па 12 
viscous terms are also of order ue d /1 ie. 9P/oy is of 
order ше а пие but the pressure difference between y пили. 
үс 4 is of order pue d /Y and the difference in дР/дх will be 


negligible compared to the external stream dynamic pressure, 


1⁄2 pue. . Therefore, for practical purposes, 
ЭР = 0 (2:10) 
ду 
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For this case, since changes in P must be of the same order 
as changes in y, the pressure does not change significantly 
through the boundary layer. 

Thus, the entire equation of motion in the 
y-direction may be dropped from further considerations. Іп 
this way, the following simplified equations are left for 


the analysis of a boundary layer: 








да + ду = 0 (225294) 
px дү 
4 
ӨРЕСІН + уби = а соли | у оц 2..1 10 
? Ox д у? д х дү 
що (22210) 
дү 


These relations are known as Prandtl's boundary layer 
equations. Unless one encounters very high Mach numbers, 
the above orders of magnitude are not changed when 


compressibility effects are considered. 


B. LAMINAR AND TURBULENT BOUNDARY LAYER 

The low viscosity fluid flow past solid bodies should be 
considered as two regions. One is the thin region near the 
boundary in which the effects of viscosity are concentrated 
and the other is the region further away from the wall in 
which the influence of viscosity is so small that it can be 
neglected. Thus, it can be stated that all viscous effects 
are concentrated in the thin region which is Known as the 


boundary layer. This boundary-layer type behavior requires 
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high Reynolds numbers. Generally, the thickness of the 
boundary layer (d) is defined as that distance from the wall 
where the velocity (u) differs by 1* from that (Ue) 
calculated by the ideal flow analysis. 

Consider a constant-property, steady, two-dimensional 
flow past a flat plate. If u/Ue is plotted against а 
dimensionless y-coordinate, 7 = (падение y 1 the velocity 
profiles are geometrically similar and reduce to a single 
curve for a laminar boundary layer flow. This is well known 
as the Blasius profile. The geometrical similarity is 
maintained regardless of the Reynolds number of the flow or 
of the local skin friction. For a turbulent boundary layer 
flow, since the viscous-dependent part and the remaining 
Reynolds-stress-dependent part of the profile require 
different length scaling parameters, there is no choice of 
dimensionless y-coordinate that leads to the collapse of the 
complete set of velocity profiles into a single curve. 

The conspicuous difference in profile shape between 
laminar and turbulent shear layers can be found in the wall 
flows. Because of the constraint on eddy size in wall 
flows, the efficiency of turbulent momentum transfer 
decreases rapidly near the wall. But, the efficiency of 
viscous momentum transfer is not dependent on y distance in 


the flow which has no heat transfer. 
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C. TURBULENCE MODELS 

The unsteady continuity and Navier-Stokes Equations are 
ЕГІС іп Есеп аш пат and türDulent flows. But, it is too 
difficult to deal with the instantaneous properties in 
turbulent flow because the turbulent flow has a complex 
time-dependent behavior. Thus, the following turbulence 
models are used to make the analysis of turbulent flow more 
convenient. 

IM P5Pandcl's Mrxting-Lengtb 

Consider two adjacent stream layers of fluid which 
move with different velocities. If a particle of fluid 
moves from one layer to the other, a momentum change occurs 
between two layers. The fast particles which enter the slow 
layer make it faster and the slow particles which enter the 
faster layer impose a drag on it. 

The mean velocity of a stream layer is u, and that 
of the other is u + 1 O9 u/gy where 1 is the distance between 
two layers. Also, the fluctuating velocity in the 
мн - кесЕјоп 15 а", and that in the y-direction is v*. 
Prandtl assumed that the turbulent fluctuations are due to 
the difference in the mean velocities of the two layers. So 
= l 3ü/əy ie. the fluctuating velocity in the x-direction 
is of the order of the difference in the mean velocities of 
two layers which have a distance 1, where 1 is the mixing 


length. Prandtl also assumed that all components of 


D 


fluctuating velocity at a given point are Of thewsameuoraam 
of magnitude. Thus, v' can be defined as v' = kl aU/ay 
where k is a constant. 

The turbulent shear stress due to momentum exchange 
between two layers is the rate of momentum transfer per unit 
area. Then the mean turbulent shear stress on the fluid is 

T= -pu'v' where u'v' %[1 2u/2y| кі да/дү| . Since the 
values of 1 and k are unknown, combine these two unknowns. 
Then, (eius о, where l is called the mixing 
length. 

2. is Cebeci оиса Model 

With the eddy viscosity concept, the momentum 

equation for 2-dimensional laminar and turbulent boundary 


layers can be written as: 





(b£")' * mtl ff" + "|1 = (Е | ЗЕ де' – 9 
2 д х 
where x is the transformed x-variable 
£ £ m/ y 


b 


2k Да 
(1 + L) (1 + Е) 


! 


Let the turbulent boundary layer бе а сопров е 
layer consisting of inner and outer regions. Then, the 


eddy-viscosity formula for the inner region is: 


2 k 
(бт); = 1 Е 
Vo 


ди 


Y — = ( 








where 1 іс the mixing length for 2-dimensional flow 
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Y 1 T CH 1 Е 
== ын exp ЕД X ( X E J 
tr | 1200 t^ 


а 
is an intermittency factor for a flat plate. 


R = ue X,p/) is the transition Reynolds number. 


Xtr 


The eddy-viscosity formula for the outer region is: 
а) 
(En), = “|| ЕЕ єр iy < y 20 
0 


Шеге 5: 2 20007 the universal constant « = 0.0168. 





Ко < 5000, е varies with Rg according to the 
the empirical formula 


< = 0.0168 125575 
TARA 


и 
A = 0.55 [1 - exp (0.243 h ^- 0.298 М) | 
h 


But this model is not used in Cebeci's interactive computer 


Rg /425 - 1 


program which will be presented in Chapter IV. 


D. TRANSITION 

The boundary layer with a finite thickness starts out as 
laminar at first in the flow past an airfoil. However, the 
boundary layer becomes unstable and all small disturbances 
begin transition to the erratically unsteady condition which 
is known as turbulence. 

In the boundary layer on blunt bodies, transition makes 
the point of separation move downstream which decreases the 
width o£ the wake. There is an abrupt change in the drag 
curve of a sphere. This change is due to a boundary layer 


effect and is also one of the transition phenomena. 


ТЭ 


The process of transition on a flatoplotem M 
incidence shows a sudden increase in the boundary layer 
thickness. Furthermore, transition involves а greatechs 1 
in the shape of the velocity distribution and vases 
decrease in the ratio of the displacement thickness to the 
momentum thickness. Also, it causes a large change in the 
skin femen. 

As described in the above, the transition from laminar 
to turbulent flow plays a very important role. Since the 
transition is not an instantaneous process and the flow is 
intermittently laminar and turbulent over a certain length 
of the airfoil, a certain point where transition takes place 
cannot be assigned. At present, there are no exact methods 
to calculate the flow within the transitional region. 
Nowadays, however, two methods - Michel's method and the е? 
method - are used to predict the transition empirically. 

1. Michel's Method 

Michel investigated many kinds of data for 
incompressible flows without heat transfer and found this 
empirical correlation, 


0. 46 
БЕ 12171 1 + 22401 | PM 


x (r 
where R9 = Је 07, в Ue x ли 
Since the momentum thickness grows more rapidly when 


the pressure gradient is positive, Michel's equation 


involves the effect of pressure gradient. However, the 
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effect of surface roughness, which is also very important, 
ШЕ ПОС included. 


? Method 


2. The e 
This е” method is o£ spatial amplification theory 
Dased on the solution of the Orr-Sommerfeld equation which 


forms the point of departure for the stability theory of 


laminar flows. 





XA а бг 
Шеге Х - х/ /”, «=, у“, В; = в; а 
ах у, S а суүріса listurbance function. 
х = 21 (Ais the wavelength of the disturbance). 
Sr is 23 сігептар frequency сш the partial oscillation. 


О те Ене amplification factor which determines the 
degree of amplification 


The basic assumption is that transition begins at 


the point where a small disturbance introduced at the 


critical Reynolds number is amplified by a factor of е” 
E. LAMINAR BOUNDARY LAYER CALCULATIONS 
I Similar solutions 
a. Blasius Solution For a Flat Plate 
Assume - a flat plate at zero angle of attack 


- steady, 2-D, laminar, incompressible 
flow. 
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- constant МмІзсевтеу 
- negligible body force 
Since the pressure along the plate is constant, 


there is no pressure gradient. The simplified Navier-Stokes 


Equation is: 





и да + у да E (2.1 
Әх дү д7 
(8706) үс 0; и = 0, у = 0 
ү =@; п = 2 


To transform from the partial differential 
equation to an ordinary differential equation, define the 


following transformation parameter: 


| 


ZU EET where += £(x,y) 
Vx 


Ф = ГЕТА Е where £ = £(1) only 


The stream function was defined to satisfy the 


continuity equation. 














ив 984 =|j/xU, ағ 07 - Ё', wheres mere (2.13) 
ду 47 ду 97 
ou c an eme UM 0 Ë 1 
ax а д y 2x ¥2 
= - Пе" (2. 18% 
2X 4 
-A 
Vus uq cy dE 27 – ЈУ Џо ИЕ 
x d7 2х 
- 1/2]YUe (£'*? -Е) (2.15) 
X 
ди = Ор df' 27 ве m (2.16 
ay а7 дү px 
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2 2 
gu = US Ue _ Шог 02 c ENS. в" (22:17) 
2У8 Ух ай ду И х 


|Ен се зө (2013) “ (2 ЦУ) into Eq. (2.12), then 


ги 


ин 
e 
` 
rh 
li 
о 


ІР су 7 ; Е! 
27 ; 


0 

о; 2 
The solution of this Blasius equation is presented in 

Msamundary Layer Theory" by Schlichting. From the 


рРапсвтогптпасіоп relation, 


x 
d = 5.0 = 5.0 d “= 177208 
X MU Л. , X [Ry 
| \/ | = 0.8604 8 = 0.664 
а (0 [RY x bs 
2 022202 ү» 
yx 
ши! = 0.664 





b.  Falkner-Skan Method 
The Falkner-skan transformation is for 2-D, 
axisymmetric laminar flow, The simplified Navier-Stokes 
Equation is: 


2 





Шош бо =; =] ЭР + y да 
дх дү fe x 9 ү? 
МЕ С) = 0. Q = 0 gos) 


у =>; u = U(x) 


Take the same n as Blasius’ büt different with a function 


Е + £(x,%) and follow the same procedure as Blasius' using 
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- 1 dP = U_dU 
P dx dx, then the transtzormeofegqu a r d T 
Ett pem T S F S ту E = сезу | = x e ЕМ = РЕСЕ 
2 д х д х 


where m is a dimensionless pressure gradient parameter 


defined by: 
m- x dU 
U dx 
(B.C) 7-0; £'-20, £ = 0 
о Е 1] 


For similarity in 2-D, laminar flow, assume f is a function 


SEE 1 only.. Then, 


Е. mè lob - ЈЕ = dnd = 0 (2.18 
2 
(В.С) 7 - 0: {f' = constant, TENE 
7 то; Е’ =1 
The fact that m is a constant leads to: 
U = с 27 


where c is also a constant. 

In the case ot m s40, 1e. U is а Е, 
Eq. (2.18) reduces to the Blasius equation. If m = 1 which 
means a 2-D stagnation flow, Eq. (2.18) becomes the Hiemenz 
equation, 

£"' « ££" - (£') * 1-20 

Some solutions of the Falkner-Skan equation for various 
values of m are presented by Cebeci and Bradshaw. 


Д Integral Methods 


a. Integral Momentum Equation 
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For steady, 2-D and incompressible flow, 





3 
т Ж = - 1 др + y du C27.) 
д х. дү ё К ay 
ou + у = 0 (029291) 
д х дү 
Пош Ед. (2.20), u ( u + day = 0 
дх дү 
At the edge of the boundary layer, О (х) 2U(x) = - 1 əP 
ð х A oX 


Then, Eg. (2.19) becomes: 


шин ee (uy) = DU(x) d U(x) + a'u 
Dx ду ах 2 y? 





Integrate this equation with respect to y, from y = 0 to 














y = , using Т = g 2u 
ду 
(M f 4 (m 
| ду = J 0(х) да дұ- | Шац оа у = - — 
‚ ах 5 2 x 0 dx PET 
* d 
Also we know, 4 = | ( = u ) dy 
0 U(x) 
0 - |" u t - u Jor (22122) 
U(x) U(x) 


Substitute Eq. (2.22) into Eq. (2.21), then, 


dU (x)0 + Г 0х) а0 (х) = Tw 
ax ax f 


This equation is known as the momentum-integral equation of 

boundary layer theory, or as von Karman's integral equation. 
b. Pohlhausen's Method 

Assume a polynomial of the fourth degree for the 


veloclty function, 
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2 З 4 
и = ЕЛ) = ал + b À + GA F ал 
U(x) 


where A is the dimensionless distance from the wall, 


оао Ди m 





(В.С) А -0; ЕЗО, 
A=1; £21, £f=0, £"-0 
Also U(x) d U(x) *y ?u -0 when Ад = 0 
ax дү” 
= 
d 9 0х) = -25 = А 
y ах 


where A is the shape factor. 
соеп F(A) = F(A) + AG(A) 


4 
2A- д А 
1/6 A (1 8 


where F(A) 
G(A) 


ж 
Thus the boundary layer parameters d, 60 and Tw can be 


determined, if the velocity profile is known, as follows: 


ж а, Ж. 
4 10 120 
2 =.1_ (92. еы 2 АРЫ 
4 63 | 5 15 144 
Ty = ШК) | 2 

d 6 


Со Thwaites' Method 


The integral momentum equation can be written 





= оу: 
а + 0 (Н +2) aUl) = е (2.23) 
dx U(x) dx 2 
(B.C) у = 0; 2u - - П(ХУК wes uo ШӘ LOK) 
g V = 3> ду 9 
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D 
Cy = Tw = » ЈЕ 
2 шим) 0 U(x) 
ЕЕ 
K = 0 40(х) 
y dx 
then Eq. (2.23) can be reduced again. 
2. 
Utx) d@ = 2 | ык) M | H(K) + 2] | = F(K) (2.24) 
y ax 


Thwaites writes an expression for F(K) 





F(K) = 0.45 - 6K 
Then we can write Eq. (2.24) as 
2 X + 
ОО 0.45 | U'(x) dx 
y 05(х) 0 


If @ is calculated for a given external-velocity 
distribution, H and С; can be determined with the following 
relations. 

КОГО СК с 0,1 


L O о КО = 1.8K 


Ве ИК — 5. 24K 


For БОЕ оо 
L = 0.22 + 1.402K + 0.018K 
0.107 +K 
Di 0.0731 + 2.088 
0.14 + K 
3. Finite-Difference Methods 


The finite difference methods are the most flexible, 
practical and efficient methods to solve the boundary layer 


equations. 
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The box method, presented by Keller and Cebeci, is 
introduced here as the preferred finite-difference method. 
The momentum equation achieved by the Falkner-Skan 
transformation can be rewritten in terms of a first-order 
system of PDE's. Then the resulting non-linear system is 
linearized by Newton's method. Finally, the block 
elimination method is used to solve the linearized 
difference equations of the boundary layer problem. 

a. Box Method. 

Using the new coordinate system, £ and n, the 
transformed momentum equation for steady, 2-D, 


incompressible flow becomes: 








£' = u (2.25а) 

ul =v (2.25b) 

(bv)! * m * 1 Ем + ти 1-и) =$fu За o 
2 3 £ Э 

(2 225 ей 





У Ы | 
Зэ 5 57 


Figure 2.1 Net Rectangle For Difference Approximation 


where £ = x and b = 1 + Vt/ y 
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A prime denotes differentiation with respect to 7. 


The boundary conditions are 


ВО 21211110) = 0 (2.26) 
Z = АЕ окш) = 1 
We denote the net points shown in Figure 2.1 by 
n n-i 
В ЕТ 2,808, n 
DEC а) 3j = 1, 2,..., J 22727) 


7 бо 
Here n and j indicate sequence numbers. 
š п 7 5 š 
With g; = gCÉn, j) denoting the value of any 
quantity at the mesh point ($n, 735), centered quantities 


can be written as two-point averages: 


п- а п 1-1 

£ -1/26$ +8 3, СЕЕ 1/20 73 «74-17 (2.28а) 
А n ща п П П 

4; - 1/2008; * $4» Ги 1/2004 * 1) (2.28ь) 


The finite difference approximations of Eqs. (2. 25а) 


ИО 2.255) for the mid point of the segment P , P are: 


2 
n n n 
Но и (2. 29а) 
h^ dia 
4 
n 
С 21 (2.29b) 
h. за 
4 


Eq. (2.25c) can be approximated similarly Бу 
centering about the mid point of the rectangle P1,P2,P3,P4. 


ЕЮ его О Eg. 12.2560) about the point 
n- A 


(5.17) 
Let the left-hand side of Eq. (2.25c) Бе L and use 
Bee (2.28b) 
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n n-i п- 2 -1 27 Я 
1/2 a s j Ë (о Е аи J - Це e ]| 
n 








K K 
VE n-i v n-! n n-i N i 
218 (ut) - (à) — (fu) Е v. — БООСОН | 
2K, 
By rearranging, 
(bv') Í m + 1 (Ем) + m 1 = сим 
2 
a- п-) 
2 Е n - 
Зу - Q0" - c£ + £” ут у" g' « (£y) | 
(bv!) * eQCEv) - u^)! « auth" Logica ан 
n- А n 
where © = 5 P та МРТ + Ay A= m^ + с 
Ky 2 
n~] -! > Е = 
Бе И на [сву ше ТЫ! 23 
L” = | (by) § + m+ 1 (fo) шэг 
2 


P 2. Centering Eq. (2.30) ароцЕе Еле ров: 
( $ 2 Де Using Eq. (2.2908 


n 
Itb]. - b; у; - SAL 
#- х h- 





2 
Тһе Еа. (2.30) becomes 
2! n n n n n 2 n 
1; (5, Va = КА Уан ) teu EVO, = ж (ш J 
jo ce пер = 2.31) 
T шин -Х д 1-Х УЕ = - ( . 
where 
new) m ` 
к и = - Па -Х + с (Ем) р-у s a 4 :3 
n-i А : * 
ШРИ = | 5; (bi yg о mti 222757 mn ж 
The boundary conditions at 3 =. аге 
EG or шэн. 
uU (2:39 
J 7 | 
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Б, Newton's Method 


For simplicity, let 23 2 а; Е 
n 
(Ер, мр, vy) at 


ВЕН оо ооу апа 2 31) сап Бе written as 


эг Гоо 24 22010122 = О 224) 
527 "аша а 2... =O (277955) 
» (bg vg - bg-, Уфа ) + (ЕГИ — 40 (ц? 4-м 
n-i | n-l n-i 
шог: РҮЭ” = Ё ¿- 1⁄4 Уа И ) = К2-ч га С) 


п-! 
Here the unknowns are on the left-hand side and Каси, 


involves only known quantities. 
Now, Newton's Method is applied to turn 
Eq.(2.33) into a linear system complemented by boundary 


Editions Eq.(2.32) and initial values 


(е) (е) (9 Р 

Жэ 0 DNE 8 боке и 
(> _ mos (е) _ "=! (2) _ 1-і 

Е) - Іі; u; = P V4 = М; ОТЕ о Sg) 
peo _ n-! (9) 2 (0) E n-! 

Е 7 = t Us 1 Ут) V5 


The superscripts in parenthesis represent the iteration 


number as follows: 


17: 1) TE 


Е; : ug Vj | i Us 122: 
и) _ „о (bp ПС (I) (1+) D с» 
25 = Е; + J fj, 0; сие дас, DNA UE Jv; 


where Jf << Е, Jax а, ду со. 


Replace f4, uz, v4 іп Eq. (2.33) with these expressions. 
Then Eq.(2.33a) becomes: 
TE a» (1) Uu) : ü) а) а) (1) | е 
АЕ: — Е, Е = 24 Е + 218200 ü j-i + J uj = 0 


3! 


Thus, 


4%; = БҮЛ 3 му т (а, | щита (2.93498 
(1) (1) (0) 
where (г); о E ћ, 4-% 
Prom Eq. О Је о Бон 
а li? 0) a (1 > 
0; + du; шоо - fus. m Pe (0 vey £ suit | = 0 
Ju; = Jugs by ЦОНХ “ул | = r3) 25 (2.34b) 
u) (D 
where (кају = uz, ` üu; + h; 
From Eg: пее) 
-i 0) n) (0) i (1) (1) 
Db TN PE Гуд) + (Е). x 2 (u 
n-i n! а; 
+ о (уру 13 27 Gee IV G- 
2 n -! ~ a сөз | у ч о) _ 201 
= Ка.) Ls [n7 (b; Va b+. MPET ) t ex, (Ту). об, (1 n 
ТЕ n =! 11) 
+ | 2 
© (узи фк 1] | | 
а ü) В (1) (1) (D) Ур a» 
here d COL, = 1/2 [tj Фу) + 2 ГЕ. fo. Дун: + ЭЙР 1 
(1) а (1) 1 0 
fu") 7-3, = » Z о; + a? du 2 
P Ry - 1/2 (d £g td £s, ) 
[un - 2 о КИ ) 
Ећеп (5,3; 49, + (ад а + (53)/ 4 6 + (5а), “Ба, + (Веи 
24 qd 4 2 4 2 д-! А 
+ is, 2 97 шу (ri), E 
- л! = Ту "m = ti) и 
where (у), = КӨРІ | hj (b? v n 27 + aC Ev) gy 
(1) һ-) (1) (1) 
-а цар тоо НЫ ғ | 
тт) оз (+) l 
` = 7 . оќ, ` = оќ 
(5,2; h; b, + оі f ex Е и. 
2 2 
-Í сі) (1) пе! 
(5,3: = - h; b, + @ f. , - < 1) 
2 2 2 4 2 2 
Uc dior (D SE: 
(5), 5 ЈЕ + = бы. 
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(21:22 - соғу. тог, 
4 — - — 220) 
4 5 2-1 2 4- 
(56), x У 233) 
0 
(p U pus 


the boundary conditions are: 


Де = 0 fu, = 0 Гај = 0 (2235) 


c. Block Elimination Method 
Eqs.(2.34) and (2.35) are the linearized 
difference equations of the momentum equation for external 
flows which has a block tridiagonal structure. We can solve 
these equations by means of the Block Elimination Method as 
discussed by Keller (1974). 
аео define themit cee dimensional vectors, 21 


and r», to express the system in matrix-vector form, 


9 


11: 


ДА 
0 Ure )3 
= |0 0236 


9 


anG “the. 3373) Mater lees: а; В 


1 0 0 го 
шэг: А; = 05 "ој mwa 
0 -1 -h,7/2 0 -1 -һ)м/2 


1 -11:122 
As = па 2 j = J 
0 
(2.32 
1 -h37/2 
Bou seg ма па 1 « Jp 
0 
0 0 0 
C4 7|0 0 0 и 
0 1 пе 
Then Eqs. (2.34) and (2.35) can be written as 
в - R (2.38) 
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where 


Ao Со do 5 
B, А, Gr Ji г, 
2 By Ag Cj 4-|4 КА 
си за. y га 
БА ~ ја b КТ 
Ву Ау L LI 
Let us factorize the matrix 6 to solve Eq. (2.38) 
(8 = хү ща 139) 
where 
I у, С. 
xr Y, C 
(5 Xj I * Y; C4 
x 4,1 23311 
х, I Y 
Here I is the 3 x 3 identity matrix. 
Xj and Уу are also 3 x 3 matrices. 
Ewecording to Eq.(2.39), we can find 
у. = Ао (2.40а) 
Жа Uc. 2, p J) (2.40b) 
PES noc Хх; су, (Ec ид JS (224 6) 
and the matrix ху has the same structure as that of the 


matrix Bç. Therefore, X4 has the elements like this, 


Фе ) 3 (Xj2)3 (х,у) 


X —- = Са). (Xa Ja (х3) 
2 | 0 4 0 2 Ü 2 
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апа y; can be denoted as: 


(Yn № (Y 00 
y e (Ya ); (Ум ЈН 
0 -1 -В 1+1 22 


From Ба: 21211007 


(ү,з о 
1 502 


1 (Уна | 
0 Cy ee) 


СҮ А 
2 


ШІ 

e 
| || 
о © 


From Ed: (2.406) 


ТЮ 
О 


21 Ux aJ 


- 17 251 (0019.9 = 0 


| | 
ІШІ 


Froma 20002) 


nou = | (Y, 5 = Е —(X 73 )3 Үз; = Aa 


уы, а (53 DEC de = (sr) = (ay з Суза Јани се оси п! (хаз; 


= -2 (s=) n 


Then we can compute the elements of x with Eqs. (2.37) and 


СОСО ro ins и 





и | (Yigg, + ha Е; У ми: Ф. | ј 
YS 2 2 
“ 2 ~ 
a ee Ce с шы]! 
| 

(хуз) = = E )4 (уз Je т (Ху )3 (Yas ger | 
оси E seal («ау бум yt S207 Ы 

Yo hy [сда СУ. о СЗС 

2 

(X - Т po. DE (8,2; О Ez Ји на guy] | 

у ( 2 = 
(X.3 ЈА == (XI y СУЕ ма t Су j (Yu ја. са (S, )4 

ie 227 2 бус, Јао 25 [9 Суу Р 
h Қал У ES E i] 
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| 


Y = 48 Dr ы 2-1 
её у -2, Жһеп х2 «В 22721) 
ШЫП, Z. = г. 22127 
27 124 — Xg Z,., Сеа о) 
where (2, ); 
24 = (2, )у (UN j 
(23). 


4 
From Eq. (2.42) for j = 0 


Тул гэ бут и в. ПД E ELE 305 

апык L <J 4 

(2,9 - (га а EE or 4 Gg, EE Xs ly (Zadar, Бх 130023745, 

(2,2; а, РЕ Тоо y (Zag, = (хз) (23)2-) 

(2; J (т) 

om Eq. (2.41), 

yy fy = 27 (2. 43а) 
уни io са 444! СУ 1) (25455) 


Then the vectors can be calculated with Eq. (2.43). The 


maree components of , for j = 0,1,2,...J-1, are: 
4 9; = — 23+ УДА E e, 
d vj I (сун), [ ); Е У | Yin Се и 52211 J 
Yi 
2 ч, - 1 | (203, Е 4 u; (у); АЛ 
(ућ 2 
where e = (23), до + Бая Vay, 
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y," аи ЕР бу 2) о. а | 


AE n (Y, P и Р Yn ја 


and for j = J 
£ uy == 2237 
j Vy == e2 (уу: л ёу зу 
(Уз ds CY ay )y = (Уз СУ n )у 
jw. = ез - (Yi y £ vz 
(у, ) y 
where e, = (2,)у - (уд )удчз 
ез = (2,3: < (уа Одит 


These calculations are stopped when 
| а | ( £ 
where v(o) is the wall shear parameter. 


6 is a prescribed value. 


F. TURBULENT BOUNDARY LAYER CALCULATIONS 

Turbulent fluid motion is an irregular condition of Е 
in which the various quantities show a random variation with 
time and space. Therefore, turbulence is characterized by 
random and chaotic motion of £luid particle sp 

The velocity varies randomly at any point ina turbulent 
fluid. The velocity components in a three-dimensional 


turbulent. fluid are: 


а с а + ц! 
v. = V ke y! 
Ww = Ww + wt 
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where u, V, W represent the mean velocities 
ta 


DNE 1 J ша апа 50 оп. 
м Lt, p 
и’. 9’, we represent the amount of fluctuation of the 


instantaneous velocities from the mean velocities. 


Obviously, we can see that 
iz 





1 |» dt - u' = O etc. 
Бут t, +, 
Similarly, the following quantities can be defined, 
Въ — 
1 jut et = u" 2 ес. 
PR LIT 
da 
n ІЗ СІН еп Vee etc. 
Бог Тт 


Even though the flow is turbulent, the time-dependent 
momentum equations are valid. However, it is impossible to 
Solve the momentum equations for turbulent flow because the 
fluctuations are random and chaotic. Thus, Reynolds 
modified the momentum equations by introducing the mean 
value and fluctating values of the flow quantities and by 


assuming the fluctuations to be continuous functions of time 








and space. From the momentum equations, 
и + о ди + у дц + ж ду = -1 Әр ЕЕ К ош ся 
at ax оу ez ^ әх ox? o y? 225 


By multiplying the continuity equations with u 


(53 + ду + au) =o 
ð X дү д 2 


39 


Sum these two equations and rearrange, then, 

















да =- 1 9P + од |у - P" д [Уд - цу] + 52 (22 - d 
әб f 3x дх| дх ду | ду 921 92 
Here take the mean value of each term using ul = 0, 
041 + ut) = au + ut = əü 
ot at ot ot 
-1 (Б + Р') = -1 Әр 
С Ə X 72 Х 
= ----- - " THER 
yu tut) - (и tut)” =)aa -— (u* + 2 uu EE 
ðX 9х 
zou о 
Ан = _ хонон Х 
y аи зала) -Уда - uv - u'v' 
ду ду 
p du + ul) -Y3ü - uw - u'w' (2.45) 
д 2 9 2 


Substitute Eq. (2.45) into Eq. (2.44) and rearrange. Then, 








PLUS а 
(б ət д Х дү 92 
= -P + ag- 0( дц + Ju'v' + dutw! 
ый СОЯН _ u u v__ u w 
2 X г 9 х дү д 2 a 


By similar procedures for the other directions, 


























ду д2 
= - ОБ «dV “A ута! + ду” + 2н") (2.46b) 
ду д ду д 2 
p í òw +u ды +V ди +w OW ) 
at ð X ду д 2 








=- ЭР tugh -ee + ды! + ix) (2.46c) 


and from the continuity equation, 








да + ду + ди -0 (22472 
aX ay д 2 
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In laminar flow, we have 4 unknowns and 4 equations. 
But in turbulent flow, we have 10 unknowns and only 4 
equations as shown at Eqs.(2.46), (2.47). This is the 
reason why we need the turbulence models which were 


discussed in Section C. 


G. SEPARATION 

In some cases, the boundary layer thickness increases 
considerably in the downstream direction and the flow in the 
boundary layer reverses its direction. The change in 
direction causes the decelerated fluid particles to move 
outwards, which means that the boundary layer is separated 
from the wall. This phenomenon, boundary layer separation, 
is always related to the formation of vortices and to large 
energy losses in the wake of the body. 

Let us consider the simplified boundary layer equation 
in order to investigate two-dimensional separation; 


+v ðu = — 1 др + Уди 


ду P X ду“ 


ü 





du 
ÒX 


and define the separation as 
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The velocity profile in the boundary layer always has an 
inflection point in the region of decelerated flow and the 
velocity profile at separation point must have an inflection 
point. Thus, separation can occur only when the flow is 
Gecelerated. In two-dimensional separation, a bubble of 
fluid which has low velocity is always formed. This bubble 
is often unsteady and distinguished by interior streamlines 
which are closed loops or which extend to infinity. 

There are two severe problems in boundary layer 
calculations when separation occurs. 


1. The Goldstein singularity at the separation point in 
direct boundary layer calculations. 


2. Numerical problems downstream of the separation point. 
If a boundary condition prescribes the pressure 

gradient, then the boundary layer methods suffer from a 
singularity due to separation. Moreover, the singularity 
causes the numerical breakdown in direct boundary layer 
methods near the separation point. This Goldstein 
singularity can be overcome by using the displacement 
thickness or the wall shear stress, instead of the pressure 
gradient, for the boundary condition, which makes it 
possible to integrate the boundary layer equations through 
the separation point. Also, the full Navier-Stokes 


equations exhibit no singular behavior. 
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On the other hand, the boundary layer equations lead to 
a numerical instability in the region of reversed flow. The 
FLARE approximation is the most common method to overcome 
this instability. The momentum transport term 
u ju/əx is deleted where u is smaller than zero. But the 
accuracy of this approximation decreases as the region of 
reversed flow increases. Therefore, it is necessary to 
introduce an upstream influence in that region. This is 
incorporated by the downstream upstream iteration procedure 
(DUIT) which consists of a sequence of alternating up- and 


downstream sweeps with the momentum transport term u du/oex. 


III. 1МТЕКАб  СП БХ МВТ ОЮ 


А. INTRODUCTION 

If the flow remains attached and the Reynolds number is 
high, the pressure distribution and the overall lift force 
can be obtained from a potential flow solution. 
Conventional boundary layer methods provide additional 
information about the skin friction distribution and tHE 
overall drag force. However, if the flow separates, no 
information 1S available for regions downstream of the 
separation point. Therefore, the overall forces cannot be 
obtained. 

For this reason, interaction methods are introduced to 
overcome this problem. They provide a special coupling 
between the inner viscous and the outer inviscid flows, and 


can be classified as follows: 


E; Weak interaction methods 
a. Direct method 
ра Inverse method 
Cc. Semi-inverse method 
2. Strong interaction method 


The weak interactions provide only a loose coupling 
between viscous and inviscid regions, i.e. two different 


regions are treated alternately. The viscous flow solver 
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deals with the flow in the viscous region and yields the 
boundary condition of the inviscid region, while the 
inviscid flow solver deals with the flow in the inviscid 
region and yields the boundary condition of the viscous 
region. The exchange of information through the boundary 
Condition is slow, but regions such as those near separation 
and trailing edges require fast and direct coupling between 
viscous and inviscid region. While the strong interaction 
method treats displacement thickness and external velocity 
Simultaneously, the weak interaction methods process one of 
these quantities as input and the other as output. 
1. Inviscid flow methods 
a. Direct boundary conditions 
(1) Prescription of the airfoil shape 
(2) деко normal velocity at the surface 
b. Inverse boundary conditions 
Prescription of a velocity distribution for 
the unknown airfoil shape 
2.  Viscous flow methods 
d DITGcNMsundgdpoccondbtlons 
(1) No slip condition requiring zero normal 
and zero tangential velocity at the 


Surface. 


45 


(2) Prescription of the external тесе а 
i.e. u-component of velocity at the 


edge of boundary layer 


b. Inverse boundary conditions 
(1) No slip conditien 
(2) Prescription of the displacement 
thickness 
с. Boundary conditions Гог simultaneous 


interaction 

(1) № slip condit on 

(2) Prescription of a linear combination 06 
displacement thickness and external 


velocity 


B. WEAK INTERACTION METHODS 
l -Direct Method 
А direct inviscid flow solver is combined with a 
direct viscous flow solver (see Figure 3.2). The boundary 
conditions are: 
ш (X70). -= v (x; Bp) = 0 (3.2) 
ü (X, v ) 280) 
This method is terminated at the point of vanishing 
skin friction. Іп the direct scheme, the pressure is 
calculated from the inviscid region, but the displacement 


thickness is determined from the viscous region. Because of 
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EN Eobenoemsnenesndeooldstern Singularity, this method 15 
Wow ao ан о отео Wiel Strong interference effects 
between viscous and inviscid regions. 
2. Inverse Method 
This method consists of an inverse inviscid and an 
inverse viscous flow solver (see Figure 3.3). The boundary 
conditions of the inverse boundary method are: 


110) (ње У(х,0) = 5. 2) 


0 
21-02) noc M Eoo] 


where 4 is WERS er eami unction. 


The roles of displacement thickness and external 
velocity distribution are exchanged in this method. The 
troubles related to the Goldstein singularity can be 
overcome, but the whoie procedure takes a long time due to 
very slow convergence. Thus, this method can be applied to 
the regions of separated flow only and needs severe 
under-relaxation. 

3. Semi-inverse Method 

A direct inviscid flow solver is combined with an 
inverse viscous flow solver (see Figure 3.4). The input is 
the displacement thickness and the output is the external 
velocity distribution for both solvers. The two external 
velocity distributions are combined and then an updated 
displacement thickness is obtained through a relaxation 


procedure. A formula for satisfactory convergence is: 
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Ж 7 
d new AS dolg (x) E B SM Ue v (z) a (3.3) 
UeI (x) 


where w is the relaxation parameter. 
The numerical weaknesses can be overcome, but the 


coupling is still loose. 


C. STRONG INTERACTION METHOD 
The simultaneous method solves the boundary layer 
equations subject to an interaction law (see Figure 3.5). 
Viscous displacement effects are allowed to cause 
substantial changes in the external velocity distribution. 
Since both displacement thickness and external velocity are 
treated as unknowns, one more additional relation, the so 
called interaction law, is needed. Thus the boundary 
conditions are: 
ПГО 050 у{х,О) = 0 
о м. P КЕ (5% 
US Cx) + 1/7 (а азе /19) 43 
X 


ч (х,у, ) 


| dí x ч 
The last equation represents the interaction law. 
Here, the external velocity for the boundary conditions 
can be written as: 
U (x,ye) = Ue (x) = Up (x) + Абе (x) (3.5) 
where Ue (x) is due to inviscid flows past the ате соли 


AU, (x) is the perturbation velocityeaue тощеше 


displacement effect of a boundary layer. 
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шэг 216 2 1 е blowing velocity concept is used. Let 
us consider an airfoil along which sources are distributed 
(see Figure 3.1). This surface distribution of sources 
represents the effect of boundary layers. There are two 
ways to predict the displacement effect of a boundary layer 


On the outer inviscid flow: 


1. Compute the inviscid flow past a displacement 
body 
2. Replace the condition of zero normal velocity on the 


surface with a condition of prescribed blowing 
velocity at the surface 


Y 


я 
V (x,d ) 


1 
| 
| 
| 
| 
| 
| 
| 





Figure 3.1 Blowing Velocity Concept 


The streamlines are displaced away from the surface by 
the distributed sources which eject the fluid at the 
surface. Theq the virtual displacement body becomes a 


streamline and the flow tangency condition is: 


x ж 
ги! (3.6) 
0.(х) ах 





From the thin airfoil approximation, the displacement 
thickness is assumed to be so small that u-component of 


velocity do not change across the layer and the airfoil in 


this connection сап be represented by а знан | - тен 
Therefore, the blowing velocity УС о) 15 еапа TE 


the source strength 6 (x). 


о с а (3.24 
Ж 
5, 4 
и(х,о) = v(x,d ) = | ду dy 
о ду 
=e” m * aUe 
- Ue dé. 2 99 
ах ах 
x 
а шулы (3.8) 
ax 
x 
Thus, а (5/7) = 172 КЕ (3. 9 


ах 
Now, the perturbation velocity due to the displacement 


effect can be written as: 


X 

AUe(x) = 1 | 04) а: (3.1% 
2s И x c 

| 


Finally, Eq. (3.6) becomes: 
о X 


Ue(x) = Ue (x) [eR dí (3.11) 
2 a 28 

Eq. (3.11) represents the interaction law in usable form and 

the integral on the right hand side is known as the Hilbert 

integral. 

In the simultaneous method, the inviscid flow solver 
provides an initial external velocity distribution, but the 
inviscid flow solver is not incorporated in the overall 
iteration process. Thus, the viscous flow solver needs the 
interaction law to obtain the necessary information about 


the inviscid region. 
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(Direct Boundary Conditions) 
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Output: External velocity 
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VISCOUS FLOW SOLVER 


(Direct Boundary Conditions) 

Input: External velocity 
distribucion 

Output: Displacement thickness 
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Comparison of old and 
new displacement 
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Figure 2.2 Direct Method 
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VISCOUS FLOW SOLVER 


(Inverse Boundary Conditions) 
Input: Displacement thickness 
distribution 
Output: Externa ооёЁ 1. 
Gis throu stem 


INVISCID FLOW SOLVER 


(Inverse Boundary Conditions) 

Input: External velocity 
distribution 

Output: Displacement body 
shape 






CONVERGENCE? 


Comparison of old and 
new displacement 
thickness 


Figure 575 Inverse Method 
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Figure 3.4 Semi-Inverse Method 
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Figuren 5 Simultaneous Method 
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то оо Кет OF^CCEBECTUSCINTERACTIVE PROGRAM 


А. INPUT DESCRIPTION 

D Intneduction 

The Reynolds number is the most powerful parameter 

which can affect the flow. Wholly different flows can 
result from different Reynolds numbers (e.g. 6.0E + 06 
against 0.28E + 06). Even though a flow with high Reynolds 
number behaves well and does not separate, it may have an 
extensive flow separation at a low Reynolds number. 
Numerical breakdown in boundary layer computations is partly 
due to this extensive flow separation at low Reynolds 
numbers. To avoid unrealistically large regions of flow 
separation, it is necessary to make some changes in the 
turbulence model for transitional flow. Reduction in 
transitional region to a proper level is a method to 
decrease numerical problems in computation. By doing this, 
more reasonable results can be obtained in low Reynolds 
number flows. Also, the computation at low Reynolds numbers 
needs more iterations to converge than at high Reynolds 
number, and the lift coefficient increases with Reynolds 
number at the same angle of attack because low Reynolds 
number flows exhibit a stronger displacement effect than 


high Reynolds number flows. 
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The results may be affected by transition location. 
The transition location can be either fixed by the user or 
computed from the following empirical formula given in 
Cebeci and Bradshaw (1977). 


0,46 
в = 1 НЫ (1 + 22400 ) Rx (4.1) 





9 = 


Laminar flow separation, however, may occur upstream 
of the transition location predicted by the above equation. 
In this case, it is assumed that the onset of transition 
corresponds to the point of laminar separation, and the 
following message will be issued in the output: 


TRANSITION LOCATION HAS BEEN RECOMPUTED AT THE 
POINT OF LAMINAR SEPARATION 


The results of this program agree well with the 
experimental results up to stall at high Reynolds number 
flows. However, low Reynolds number flows may not agree 
well or experience numerical breakdown close to stall with 
the following messages: 

at the very top of output, 

+ LEY 002 2 Si ORs 
or at the very bottom of output, 


IFY 207 I VFNTH: PROGRAM INTERRUPT (2)- 
FLOATING-POINT EXCEPTION OVERFLOW 


IFY 259 I SINCT: /ARG/ = Zar güne ns ЕЕ 
(HEX = hexadecimal), APPROACHES SINGULARITY 
or MAXIMUM NUMBER OF ITERATIONS EXCEEDED 


Thus, the range of angle of attack should be considered 


carefully. 
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If the user requires many sweeps with IPRNT = 0 or 
2, then the user needs to choose class "G" or "J" for enough 
running time and add additional command ", LINES = (m)" just 
after nnnnP on the following line to make enough space for 


Print: 


11% MAIN ORG = МРСУМІ nnnnP INES - (mJ 


where nnnn is user's ID number. 
m is the number of output lines in thousand 
(Generally 10 ys enough) 
If the user does not do that, the following error 
message will be issued at the very top of the output and 
printing will be stopped abnormally: 


IAT 1600 JOB number (user's job number) LINES EXCEEDED 
IEF 450I user's job number - ABEND S722 UO0000 


2. Detailed Input Data Description 
The input to the computer program consists of 1 
title line, 3 control lines and airfoil coordinate. User 


must follow the data format for specified column and type. 


Title Line | FORMAT (18 A4) 


This line provides any description as desired with 


any acceptable machine characters. 


ог го Line 1 | FORMAT (5 I5) 


— w w x w x  —  — — w  — n инээ 


2) 


Column 


5 


10 


185 


20 


25 


Name 
ITE T) 


ЇЇВӨ2) 


IRSTRI 


Read Starting Solution 


No 
Yes 


Explanat rem 


Transition flag for the 

lower surface 

Transition flag for the 

upper surface 

= 1 Point of transition 
has to be specified 
by user. This opto 
shall be used if ex- 
perimental results 
are available or to 
avoid osciilations of 
the computed transition 
points. In the latter 
case, fix transition at 
the most upstream 
occurrence of computed 
transition. 

= 4 Point of transition 
will be calculated 
according to Eq. (APEE 
No input is necessary 
for XTRL and ХТЕО. DES 
computedpoint o£ 
transition will be 
redefined, if laminar 
separation takes place 
upstream of the 
location predicted by 
Eq TS 
Boundary layer restart 
flag. 


Store Final 
solution or Unite 


Yes 
Yes 


Leave 0 for general use 


ТСЬМАХ 


ІРЕМТ 
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Number of sweeps. 


Low Reynolds number 
flows need more sweeps 
Lo converger 


Print flag to cont rom 
output printi 


= —1 Summary print for the 
last two sweeps 


= 0 Summary print for each 
sweep. 


-277 Whole print for all 
Sweeps. 


w control Line 2 | FORMAT (4 E10.0) 


— — — => > => = > с о  — — — 


Column Name Explanation 
0 RL Chord Reynolds number 


uot ATRO Fixed transition location 
for lower surface 


21-30 ХТКО Fixed transition location 
for upper surface when the 
stagnation point is above 
the leading edge, XTRL may 
be positive. When the 
stagnation point is beneath 
the leading edge, XTRU may 
be negative. 


91-40 ALPO Angle of attack in degree 
Mmeeontrol Line 3 | FORMAT (15) 
Column Name Explanation 
Nm MPTS Number of airfoil 
coordinates. 
mecoordinate Data Line | FORMAT (2F 10.0) 


Input dimensionless airfoil coordinates as two columns in 
опе row ( x/c, y/c). The order is trailing edge — lower 
surface — leading edge — P upper surface —— trailing 
edge. Thus, trailing edge is input twice. 
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B. OUTPUT DESCRIPTION 
1. Introduction 

To take the appropriate results from the output, 
user must check the convergence. This can be done by 
comparing the convergence indicators, lift coefficient and 
displacement thickness at trailing edge, when the 
computation is completed successfully for the given sweeps. 
If the convergence indicators show steady values over some 
sweeps i.e. each difference of the convergence indicators is 
less than 1%, the results are considered as converged. 

Sometimes, the user may experience failure in 
computations at low Reynolds numbers and high angles of 
attack due to numerical breakdowns. These breakdowns take 
place when Newton iteration does not converge near to the 
trailing edge on the upper surface or the computation is 
terminated by Fortran error with this message: 

IFY 2511 SSQRT; ARG = argument, LESS THAN ZERO 

Also, extensive flow separation cause numerical 
breakdowns of the boundary layer computation. These 
unrealistically large regions of separation at low Reynolds 
numbers can be reduced by changing the transitional flow 
model. 


27 Detailed Output Data Description 


The output of the computer program for IPRNT - -1 
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and 0 can be divided into three parts as follows: 
КИ та ала ии 1 ела lift coefficient 


2. Alternating boundary layer parameters for each 
surface and wake. 


3.  Inviscid and viscous pressure distributions 
The following is the list of computed boundary layer 


parameters in the order in which they appear in the output 




















meint. 
Name Name 
ШИП рг тп їп Drogram Definition Explanation 
ACNX) 2441) r Dimensionless sur- 
С face distance from 
the stagnation point 
X/C(NX) ХЕСТ) X Dimensionless chord- 
С wise distance from 
leading edge 
У(1,МХ) же (42%) Dimensionless wall 
нг shear stress para- 
meter. 
ЕЕ СЕ(І) 2 Tw Босат етш friction 
P Ue coefficient 
X 
DELST DLS(I) Z Dimensionless dis- 
С placement thickness. 
* 00 
Я = | (1 T PPS | dy 
о Џе 
ТНЕТА ТНЕТА (1) Oi Dimensionless 
ра C momentum thickness 
6 4 ü (1- ü ы 
T Ue Ue 
UECNX) UE(CI) UeI Dimensionless exter- 
О nal velocity 


computed from the 
inviscid method with 
viscous effect. 


W(NP,NX) WNP(CI) Ue V 


D(NX) D( KX) x 
Uev d (вт. 
Ü C 

DB(NX) DB( KX) 


а ШОЕ А (мх) | | + @ Í 


IT ITN(I) 


NP NNP(I) 


C. HOW TO CHANGE THE ORIGINAL PROGRAM 


Dimensionless exter- 
nal velocity 
computed from the 
interactive boundary 
layer method. 


Product of veloce ин 
and displacement 
thickness in dimen- 
sionless form for 
the current sweep 
(D), £rom the pre- 
vious sweep (DB). 

D and DB are printed 
in the output aspa 
dated values through 
the relaxation 
procedure 


W (NP, Nx) |) | 


UE (vx) 


Iteration count for 
the convergence of 
solutions at a given 
NX-station. 


Number of points 
across the boundary 
layer. 


Cebeci's computer program is written entirely in FORTRAN 


and consists of the following eight FORTRAN files: 


FILE 1 FORTRAN Al 


EICE 2 FORTRAN Al 


FILE 3 FORTRAN Al 


FILE 4 FORTRAN Al 


ThEE———FORTRAN—RI 


EILE 6 FORTRAN Al 
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ЕЕ И 


ШЕК ӘЗ 


FILE 2 


FORTRAN Al 
FORTRAN Al 


FORTRAN Al 


Also, the user needs the following six JCL files to run 


changed program: 


OUEST 1 
ALLPDS 
QUESTS 
STOPDS 
LOADMOD 
INTAPCLG 


ШЕСЕЗТ 1" 


КЕЕЕРр5 " 


UEST 3" 


I TOPDS" 


JCL А1 
JCL А1 
200 А1 
JCL A1 
JCL Al 
ICL àl 


lists all data set on MSS (Mass Storage 
System) that belong to student user ID 
number and in group PUBAB (public disk 
volume), which deletes the data set 180 
days after creation. 

allocates space for a PDS (partitioned 
Data Set) on PUB4B prior to reading a 
data set into it, using the utility 
program IEFBR 14 which pre-allocates or 
deletes data sets. 

lists the members of a PDS and calculates 
the remaining available space within the 
data set. 


places a FORTRAN source file ina PDS. 
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The SYSUT2 DD statement describes the out- 
put data set and the SYSUTI1 DD statement 


describes the input data set. 


"LOADMOD" creates a load module. 


"INTAPCLG" coinpiles source.code, leads and executes text file. 


A list of all JCL files is given at the end of this chapter 


and more information is in the User's Guide to MVS at NPS 


C1986); 


The following is the procedure to change the program 


апа гип lt: 


Step 1 


p 


Make changes in FORTRAN files as required. Remember 
those eight FORTRAN files compose one program. 

Thus, the user must check all files which are 
related to changes. 


2. Review and update FORTRAN errors, if necessary. 

Step 2 

1. Change job name "xxxxxx", user ID number "nnnn"' апа 
library name "yyy" with user's own in all JCL files 
to make Chem be the user's. 

2. "submit four JCL files: 
"OUESTI"  "ALLPDS" "OUEST 3*5 2 1 22 

Step 3 

Do the following for all FORTRAN files to submit 

them to user's library. 

1. Open "STOPDS" JCL file. 

2. Change the FORTRAN file number indicated by n on the 


sixth line. 


//SYSUT2 DD DISP === УТЕРИ ИЕ 
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3. Insert the FORTRAN file between 
ОТ: рр 5 апа ние". 


DEM ели ел IOPDS" JCE files 


5, Check the result by any error message at the very 
ората осо Part ald condition code. 


6. Update errors, if necessary. 
Step 4 
Do the following to run the changed program. 
ша брег INTAPCLG" JCL file. 
2. Insert the input data between 
о АІ О ж апа ?/%". 
Е Submit "INTAPCLG" JCL file 


From the second change, do steps 1, 3 and 4. 


D. APPLICATION OF CEBECI'S PROGRAM 

Cebeci's interactive program was applied to a single 
airfoil, FX 63 -137, at three Reynolds numbers. The print 
flag is O and the number of airfoil coordinates is 49 for 
all cases. The results are compared with experimental data 
and the turbulence model has been changed to get better 
results. The experimental results are taken from Reference 
Br. 

Figure 4.2 and 4.3 show the comparison of lift and drag 
coefficients for three Reynolds numbers. Drag coefficients 
are computed at the trailing edge. It is seen that the 
results are closer to the experimental data as the Reynolds 


number increases. 
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Figures 4.4, 4.5 and 4.6 show the ӘКІП сеје ша 
coefficients on the upper and lower surfaces as a function 
of angle of attack for the same three Reynolds numbers. 
Legends indicate the angle of attack. As the Reynolds 
number increases, the length of the separation bubble 
decreases on both surfaces for the same angle of attack. 
Especially, Figure 4.4, which has a different label of 
y-axis from the other figures, indicates an unrealistic P i 
on the upper surface at low Reynolds number. In Figure 4.5 
and 4.6, however, the flow on the upper surface is separated 
from about 80% chord at all angles of attack. 

Figures 4.7, 4.8 and 4.9 present how the displacement 
thickness varies according to the angle of attack and the 
Reynolds number. As the Reynolds number increases, the 
displacement thickness decreases on both surfaces. As the 
angle of attack increases, it also increases on the upper 
surface, but decreases on the lower surface. 

Table 4.1 represents the computed and fixed transition 
locations which are used in Figures 4.10 4.13. The 
computed transition locations XTRL = 0.288922 at ALPD = 8 
degrees and XTRL = 0.206866 at ALPD = 10 degrees were 
printed incorrectly. In this case, the correct transition, 
locations can be obtained by the following procedure: 


l. Calculate the difference of XTRL's at ALPD - 4 
degrees and 6 degrees. 


2. Add the difference to XTRL of ALPD = 6 degrees and 
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take the result as a temporary XTRL for ALPD = 
8 degrees. 


СООК тање ком Location with the temporary XTRL 


Dune ues pogramubesprinte (IJ, GBMIRGIO from 
SUBROUTINE OUTPUT in FILE 3 FORTRAN Al. 


5. Take XC(I), where the value of GAMTR(I) is finally 
equa letomzero, as the correct transition location. 


6. From the same procedures with XTRL's at ALPD = 6 
degrees and 8 degrees, obtain the correct transition 
location for ALPD = 10 degrees. 

Figures 4.10~ 4.13 show the effect of variations in the 
empir Ical constant, 507 , in turbulence model at Re = 0.28 
X 10° . In the legend of Figures 4.10 and 4.11, the three 
numbers are the values of the empirical constant, C means 
ес the transition location is computed by Eq. (4.1) and F 
means that it is fixed. The lift curve is closer to the 
experimental data as the empirical constant decreases, but 
the drag curve is closest when the empirical constant is 
120. 

The reason why we can obtain better results by reducing 
Ene empirical constant is shown in Figure 4.12 and 4.13. 
GAMMATR (Yip) is the intermittency factor which is a 
function of the x-coordinate with values 0.0 at the 
beginning point of transition and 1.0 at the ending point of 
transition. The intermittency factor makes it possible to 
avoid a sudden transition from laminar to turbulent by 
smoothing out the step-shaped change of viscosity from 


Kinematic to eddy. The relation between the empirical 
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constant and the intermittency factor is given as: 
г = l - exp - Е - Хүр! Р Re x. | I dx | 
Yer 227 


where, Xr denotes the beginning point of transition. 





The transition length is mainly determined by the 
empirical constant y which is set as 1200 in Cebeci's 
original program. Decreasing the value of the empirical 
constant reduces the transition length. In this way, the 
unrealistic £low shown in Figure 4.4 can be avoided and more 
reasonable results, which are closer to the experimental 
data, can be obtained. In Figure 4.12 and 4.13, the legend 
indicates the angle of attack, and each of four line 
patterns is used twice for two angles of attack (e.g. — ; 

- 4 degrees and 4 degrees). Тһе direction of curves; asm 
angle of attack increases, is from right to left at the 
upper surface (Figure 4.12) and from left to right at the 


Lower surface (figure 4. 15). 


68 


Јоле Блок је Files 


| 


XXXXXXX ALLPDS JCL Al ххххххх 


//Z XXXXXX JOB CNNNN,9999), 'ALLOCATE PDS',CLASS-A 
//XMAIN ORG=NPGVM1] .NNNNP 

22 EXEC PGM=IEFBR14 

Z//SYSPRINT DD SYSOUT-A 

//DD1 DD UNIT=3330V,MSVGP=PUB4B,DISP=(NEW,CATLG, DELETE), 
17 SPACE=(CYL,(4,6,6)),DSH=MSS.SHHHNHN.YYYLIB 

/” 


XXXXXXX QUEST] JCL А1 kxxxxxxx 


722999994 JOB (NNNN,9999)5, 'QUESTION', CLASS-A 
//¥*¥MAIN ORG=NPGVM1 .NNNNP 
// EXEC PGM=IDCAMS 
//SYSPRINT DD SYSOUT=A 
//SYSIN DD X 
LISTDSET GROUPCPUBGB) LEVELCMSS.SNNHN) 
/ X 
/” 


XXxxxxx  QUEST3 JCL Al ххххххх 


// XXXXXX JOB (NNNN,9999), 'LIST', CLASS-A 
//¥*¥MAIN ORG=NPGVM1 .NNNNP 
22 ЕХЕС РОМЕТЕНЕТЪТ 
//SYSPRINT DD SYSOUT=A 
//DD1 DD UNIT=3330V,VOL=SER=MS0005,DISP=SHR 
//SYSIN DD х 
LISTVTOC FORMAT,VOL233350V-zM350005, DSNAME-MSS.SNNNN.YYYLIB 


LISTPDS VOL=3330V=MS0005, DSNAME=MSS .SNNNN.YYYLIB 
/ X 
Ж 


XXXXXXX STOPDS JCL Al kxxxxxxx 


// XXXXXX JOB (NNNN,9999),'PLACE PDS' 

//XMAIN ORG-NPGVMl.NNNNP 

27 EXEC PGM=IEBGENER 

//SYSPRINT DD SYSOUT-A 

Ze oY SIN DD DUMMY 

22515012 DD DISP=COLD,KEEP), DSNAME=MSS.SNNNN.YYYLIBCINTAIR_), 
77 DCB=(RECFM=FB,LRECL=80,BLKSIZE=8000) 

5750971 DD X 


/ X 
/” 
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xxxxxxx LOADMOD JCL Al  Xxxxxxx 
УЛХХХХАХ 

//ЖМАЈТН ORG=NPGVM1 .NNHNP 

// EXEC 


Z FORT. STSERINT DDT DUMNY 


72 КТ 20013 DD DISP-SHR,DSHSMSS. 
⁄ 7 DD DISP=SHR,DSN=MSS. 
17 DD DISP=SHR,DSN=MSS. 
/7 DD DISP=SHR, DSN=MSS 
27 DD DISP=SHR, DSN=MSS 
// DD DISP=SHR,DSN=MSS 
22 DD DISP=SHR, DSN=MSS. 
22 DD DISP=SHR,DSN=MSS 


SNNNN 
SHNHH 
УНИНН 


. SHHNH 
„ЭНИНИН 
. SNNHH 


ОННЫН 
SHANNA 


JOB (NNNN,9999), "CREATE LOADMODULE',CLASS=C 
FORTVCL,PARM.FORT=*LVLO77),NOS,NOX,NOMAP * 


. YYYLIBCINTAIRI) 
.YYYLIBCINTAIR29 
.YYYLIB(INTAIR3) 
.YYYLIBCINTAIRA) 
.YYYEIBCINUABERSE 
-YYYLIBCINDATRZ? 
A и 
9 


YYYLIB(IHTAIR 


//LKED.SYSLMOD DD DISP=SHR,DSN=MSS.SNNNN.YYYLIBCINTCAS) 


/7LKED.SYSIN DD DSH-NULLFILE 
Уу 


жжжжжжж INTAPCLG (С. Al хххххж 


/ / XXXXXX 
//XMAIN ORGZHPGVMI.NNHHP 
// EXEC 
// REGION.GO=1280K 
//FORT.SYSPRINT DD DUMMY 


7/ FORT. SYSIN DD DISP-SHR,DSN-MSS. 
27 DD DISP=SHK;DSN-HBSS. 
27 DD DISP=SHR, DSN=MSS. 
77 DD DISP=SHR, DSN=MSS. 
722 DD DISP=SHR,DSN=MSS. 
// DD DISP=SHR, DSN=MSS. 
// DD ОТ5Р=5НК,П5Н=М55. 
22 DD DISP=SHR, DSN=MSS. 
ЖЕКЕ. ЗҮЗЕКТИТ DD DUNGY 


ZZLKED.SYSIN 
//60.ЕТ01Ғ001 
ZZ60 FTOZFUOIT 
2760. РТОЗБЕООЈ 
/760.ЕТ04Ғ001 
//G0.FT99F001 
//60.ЕТ08Ғ001 
//60.ЕТ09Ғ001 
РУСО ЗЕТТОЕ ВО 
“ЖОО ЕТ 1 ТЕВО 1 
Z/G0.FTOSFODI 
//60. SYSIN 

/ Ж 

22 


DD DSN=NULLFILE 


UNIT=SYSDA,SPACE=(C 


SYSOUT=A 
х 
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UNIT=SYSDA,SPACE=(CYL, 
UNIT=SYSDA,SPACE=(CYL, 
UNIT=SYSDA,SPACE=(CYL, 
UNIT=SYSDA,SPACE=(CYL, 
UNIT=SYSDA,SPACE=(CYL, 
UNIT=SYSDA,SPACE=(CYL, 
UNIT=SYSDA,SPACE=(CYL, 


UNIT-SYSDA,SPACE-CCYL, 


x 


SNNHN 


SHHHH. 
SNNNN. 
SNNNN. 
УНННН. 


SNHNN 
SNNHN 
SHHHH 


YL, 


(1,1)) 
(1,1)) 
(1,1)) 
CI 1925 
(1,199 
(1,1) 
(1212) 
(1,1)) 
(1,1)) 


JOB CNNNN,9999), 'COMP LOAD AND GO',CLASS76 
FORTVCLG,PARM.FORT-'LVLC772),NOS,NOX,NOMAP', 


.YYYLIBCINTAIRI1) 
YYYLIBCINTAIR2) 
YYYLIBCINTAIRS) 
YYYLIBCINTAIRG) 
YYYLIBCINTATR6) 
.YYYLIBCIHTAIR?7) 
. YYYLIBCINTAIR8) 
.YYYLIBCINTAIR9) 
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V CONC EOS TONS 


Cebeci's interactive computer program was applied to 
the Wortmann-Althaus FX 63-137 airfoil to show the 
pupa vis stpongovulscous/inviscid interaction methods to 
predict airfoil flows at low Reynolds numbers. 

From the comparisons with the experimental data, it 
was confirmed that the results are closer to the 
experimental data as the Reynolds number increases. 

Also, much better results were obtained by decreasing 
the empirical constant uL. 

Therefore, it was concluded that the boundary layer 
transition model has an important influence on the 
pucuuctive Capability of viscous/inviscid interaction 


methods. 
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